
f-HNICAL REPORT SECTfOW 



NPS-69BC76101 



NAVAL 

>) 



POSTGRADUATE SCHOOL 



Monterey, California 



DIFFERENCE TREATMENT OF TRANSIENT 
IRE DISTRIBUTION IN AN INSULATED PIPE 

by 

John E. Brock 
October 1976 




FEDDOCS 

D 208.1 4/2: NPS-69BC761 01 



NAVAL POSTGRADUATE SCHOOL 
Monterey, California 



Rear Admiral Isham Linder T D n 

Superintendent t' R ’ Bors ting 

Provost 



FINITE DIFFERENCE TREATMENT OF TRANSIENT 
TEMPERATURE DISTRIBUTION IN AN INSULATED PIPE 

Wiis report discusses a Crank-Nicolson type of finite 
difference formulation for the problem of transient 
temperature distribution in a uniform cylindrical pipe 
w ^T h external insulation and containing a fluid for 
which temperature, pressure, and flow rate ere given 
as functions of time. All material properties aro 
assumed to depend on temperature. Discussion is also 
given of a digital computer program which implements 
this formulation. This implementation incorporates 
properties of superheated steam, the internal fluid, 

^7 1//4 Cr 1 Mo ( p 22), the pipe material. Stress 
calculations are also treated. 



NPS-69BC76101 
October 1976 



UNCLASSIFIED 



SECURITY CLASSIFICATION OF THIS PAGE (When Dele Entered) 



REPORT DOCUMENTATION PAGE 


READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


1. REPORT NUMBER 2. GOVT ACCESSION NO. 

NPS-69BC76101 


3. RECIPIENT’S CATALOG NUMBER 


4. TITLE (end Subtitle) 

FINITE DIFFERENCE TREATMENT OF TRANSIENT 
TEMPERATURE DISTRIBUTION IN AN INSULATED PIPE 


5. TYPE OF REPORT & PERIOD COVERED 


6. PERFORMING ORG. REPORT NUMBER 


7. AUTHORfs; 

John E. Brock 


8. CONTRACT OR GRANT NUMBER^ 


9 PERFORMING ORGANIZATION NAME AND ADDRESS 

Professor John E. Brock (Code 69Bc) 
Department of Mechanical Engineering 
Naval Postgraduate School 


10. program element, PROJECT, TASK 
AREA 6 WORK UNIT NUMBERS 


11. CONTROLLING OFFICE NAME AND ADDRESS 


12. REPORT DATE 

October 1976 


13. NUMBER OF PAGES 


14. MONITORING AGENCY NAME 6 ADDRESS^// different from Controlling Office) 


15. SECURITY CLASS, (of thta report) 

Unclassified 


15a. declassification/downgrading 

SCHEDULE 


16. DISTRIBUTION STATEMENT (of this Report) 

Approved for public release; distribution unlimited. 


17. DISTRIBUTION STATEMENT (of the abstract entered in Block 20, If different from Report) 


18. supplementary notes 




19. KEY WORDS (Continue on reverse side If necessary and Identify by block number) 

Transient temperatures in pipe 
Transient radial temperature 
Nonlinear temperature transients 


20 ABSTRACT (Continue on reverse side If necessary and Identify by block number) 

This report discusses a Crank-Nicolson type of finite difference formulation 
for the problem of transient temperature distribution in a uniform cylindrical 
pipe with external insulation and containing a fluid for which temperature, 
pressure, and flow rate are given as functions of time. All material proper- 
ties are assumed to depend on temperature. Discussion is also given of a digi- 
tal computer program which implements this formulation. This implementation 
incorporates properties of superheated steam, the internal fluid, and of 2-1/4 
Cr 1 Mo (P22), the pipe material. Stress calculations are also treated. 



DD , jan M 73 1 473 EDITION OF 1 NOV 65 IS OBSOLETE UNCLASSIFIED 

S/N 0 102-014- 6601 | 



SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 



FINITE DIFFERENCE TREATMENT OF TRANSIENT TEMPERATURE 
DISTRIBUTION IN AN INSULATED PIPE 

by J. E. Brock 

Abstract: This report discusses a Crank-Nicolson type of finite 

difference formulation for the problem of transient temperature 
distribution in a uniform cylindrical pipe with external insulation 
and containing a fluid for which temperature, pressure, and flow 
rate are given as functions of time. All material properties are 
assumed to depend on temperature. Discussion is also given of a 
digital computer program which implements this formulation. This 
implementation incorporates properties of superheated steam, the 
internal fluid, and of 2-1/4 Cr 1 Mo (P22), the pipe material. 
Stress calculations are also treated. 

CONTENTS 



Abstract 

Contents 

Body of report 

Objective 

Derivation of governing ML PDE 

Finite difference formulation 

Method of solution 

Digital computer program 

Properties of the solution 

Appendix A - Discussion of the subroutines . . . 

Appendix B - Stress calculations 

Appendix C - Instructions, data preparation, etc. 

Appendix D - Program verification 

ADpendix E - Heat transfer at outer radiug . . . 

References 

Initial Distribution List 



1 

1 



. 2 
. 2 
. 3 
. 5 
. 6 
. 6 
. 10 
. 13 
. 19 
. 24 
. 28 
. 32 
. 34 



- 1 - 



1. Objective 

The purpose of this monograph is to present and to document a 
procedure for determining transient temperature distribution and the 
corresponding significant thermal stress components in metallic pipes 
which are insulated externally and which contain steam the tempera- 
ture, pressure, and flow history of which is specified. Although a 
classical solution is known (13)* for an idealized version of this 
problem, it is not canable of dealing with the variable heat trans- 
fer coefficient which governs the exchange of heat between the oipe 
wall and the flowing steam it contains; this coefficient is a comp- 
licated function of the properties of steam. Furthermore, the an- 
alytic solution presumes constant metal thermal conductivity and 
thermal diffusivity while actual piping materials exhibit very sig- 
nificant variations of conductivity and diffusivity over ranges of 
temperature usually encountered in engineering practice. Accordingly 
a numerical approach, capable of dealing with these variations, has 
been chosen for the procedure. 

2. Derivation of the governing nonlinear partial differential equation 

A convenient point of departure for the derivation of the govern- 
ing nonlinear partial differential equation is the relation ( 5 ) 



, 9T a /Uv. 9T, . 9 ,yA. 9T* . 9 .Xu 9T, 

Xuvoo 5t = 3t ( r k 3s ) ♦ 5nW + 



(i) 



Here and later T denotes tenperature , t denotes time, c denotes spec- 



* References are listed on pages 28 and 29. 
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ific heat at constant pressure, and p denotes mass per unit volume. 
Here £, n, and z, denote generalized space coordinates and A, p, and 
v are coefficients such that the distance ds between points (f, ,n,f) 
and ( F,+d£,n+dn ,C+d£) is given by 

(ds) 2 = X 2 (d£) 2 + u 2 (dn) 2 + v 2 (drj 2 (2) 

Specializing to the cylindrical coordinate system, we take t, = r, 

JVH 

n = 6, C = z. Furthermore, we assume axial symmetry so that ~ = 0, 

and we presume no variation with respect to axial position so that 
9T 

•gr = 0. We thus easily arrive at the evaluations A = v = 1, ii = r, 



and at the equation 
1 9T _ d 7r . 



. 1 5T 13k 3T 

a 9t 9r z r 9r k 3r 3r 



(3) 



where a = k/pc is the diffusivity . The variation of k with respect 
to r is due solely to the fact that conductivity is a function of 
temperature which itself depends on radius. Thus we finally arrive 

at 1 9T _ 9 2 T . 1 9T , ,k' 97 2 

a 9t 9 r 2 r 9r k 9r 

where k* = dk/dT for the pipe material. If the variation of k with 
respect to T were simple, there would be some point in dealing with 
the Kirchhoff transformation ( b ) . However, it is not, and we will 
deal directly with the nonlinear partial differential equation 
3. Finite difference formulation 

We employ a Crank-! Ji cols on finite difference formulation, as 



follows. Equation 4 is differenced in time according to 
C(HT) + + (HT)]/2 = (T + - T)/[5(« + + «)/ 2] 



(5) 



where 6 represents the time increment, superscript ( ) indicates 
evaluation at the advanced time t + 6 , and H is the space operation 
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( 6 ) 



9 2 T . 1 9T 



91 \ 2 



W = “4 + + (k'A)(~) 



9r' 



The latter is next represented approximately by a finite difference 
operation. We use the notations 



9 2 T 



A = radius increment, A = — =- 

9r 2 



R _ 9T 

» B " 9? 



, p = k'A 



(Note the new use of symbol p.) Then, generally, 

T 1+1 = T ± + BA + AA 2 /2; T = T -BA + AA 2 /2 (7) 

where subscript (^) indicates evaluation at r = r. = n i , subscript ) 

indicates evaluation at r = r i _ 1 = (n^-1) , etc. (In general n^ will not 
be an integer.) We easily find 

A = (T i+] _ + T._ 1 -ZT^/A 2 ; B = (T i+] _ - T-^AA (8) 

so that 

«T - (T^rt^-2^ ♦ + » 1 (T 1+1 -T 1 . 1 )/‘0)/A 2 

Setting = 4A 2 /6(at+op ) , we arrange the general equation in the form 

[-1 + — + p!(T + _ l 1 -T + .)/i|]T! . + (B.+2)Tt 
2n^ l l+l l-l l-l i i 



(9) 



- [1 + sr + p i (T i+r T i-i )/, ' :iT i 



i+1 



- + T in + T t-i + (T i + r T i-i )[ S7 + p i (T i+r T i-i )/,,] (10) 

On the left, the bracketed coefficients of the advanced temperatures 
Tt_i and T+ + ^ themselves contain these (unknown) advanced temperatures. 
Accordingly, it is contemplated that iteration must be used in which earlier 

-}• -j* 

evaluations of T^_ 1 and T^ +1 are used in evaluating the bracketed coeffi- 
cients . 

Equation 10 is valid for interior nodal points, i = 2,3,...,N, where 
N is the number of subdivisions of the pipe wall. At the inside surface 
where the radius is r^ = n^A, there is a convective boundary condition 
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*i!| = WT.-*) 



(ID 



where <J> Is the fluid terrperature and h is the surface heat transfer 
coefficient. We write this as 

= B = (o/A)(T -<fr) (12) 

r i 

where o = hA/k. , and we also use the first of equations 7, taking 
i = 1. Thus we get 

HT = {2(T 2 -Tj ) - o(T 1 -<D)[2 - - o x a{T r <b) ]}/A 2 (13) 

Thus, the appropriate equation for the inner surface is 
+ 2 + o + [2 - - p^a + (Tt-2d + )]}Tt - 2T+ 

= (3,-2)T, + 2T 2 - o(T 1 -d) [2 - - o x o(T -4)1 

n 1 1 

+ a + <t> + (2 - pj— + p*o + <}> + ) (W) 

where 6 t = 1 TA 2 /6(a^+a j ) . Again, note the nonlinearity represented 
by the appearance of the advanced temperatures in the bracketed coef- 
ficients and on the right side. The advanced terrperature also appears 
in B j and in o + since a* and h depend on the advanced inside surface 
metal temperature. 

Finally, the last equation pertains to the exterior surface node 

3T 

(numbered M = N+l) where there is perfect insulation so that = 0. 

Employing also the second of equations 7 with i = M, we get 

ITT = 2(T n -T m )/A 2 

and the last difference equation becomes 

< V 2 > T S - 2T S - ( 8 r 2 )T n + 2t h 

A. Method of solution 

Because of the nonlinearities, an iterative solution is used. 
However, since the coefficient matrix is tridiagonal, it is convenient 
and expeditious to use a pivotal method several times rather than to 



(15) 

(16) 



3T 

3r 
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use a purely iterative method sucn as that of Gauss-Seidel. V/e nay 
think of the equations in the form (cf. Reference 1*1 ) 





a i T t-i + b i T t + ’ d i (1=1 - 2 -- 




(17) 


with a^ = c.^ = 0. 


Then we form 






o 

11 

*o 

o 


c i = c i^ b i“ a i c i-l^ (i=2 , . . . ,N) 




(18a) 


• ft 

o 

II 

*o 

T) 


d i = (d i" a i d i-l )/(b i" a i c i-l ) (i=l,2,. 


. . ,M) 


(lfib) 


and obtain the advanced temperatures as 






• ft 

It 


T i = d i ” ^i+l ( 1=N * N_1 » • • • » 2 • D 




(19) 



Because of the appearance of the advanced temperatures in the coef- 
ficients, several iterations are required. 

5. Digital computer program . 



The digital computer program discussed in appendix C hereto 
implements the theory developed above. The main program is called 
PIPE1EM. The properties of the fluid are determined by subroutines 
CSUBP , PRAND, and ABSMJ. These are called by subroutine AITCH which 
determines the surface heat transfer coefficient h. "Tie only fluid 
whose properties are incorporated in the present program is superheated 
steam. The properties of the metallic piping material are determined 
by subroutines METAL, COEXP, and EYUNG. All these subroutines will 
be discussed further in appendix A hereto. 

6. Properties of the solution . 

Having determined the advanced temperatures T^, (i=l,...,M), 
certain significant properties of this solution are determined in a 
subroutine called AVERT. These properties are TAV, DTI, UV 2, SL, 
and SR which will be defined and described at this point. 
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TAV, DTI, and UV2 are rigorous counterparts to "average temper- 
ature," "delta- tee- one," and "delta-tee-two" given in various piping 
and vessel codes for the estimation of thermal stress. As will be 
pointed out in appendix B, for uniform cylinders with axial symmetry 
of loading and temperature, the proper stresses to be used in code 
evaluations may be determined directly without the use of these quan- 
tities; however, their use is widespread and they are calculated for 
conparison purposes. The present writer is responsible for the defi- 
nitions given in the codes. For simplicity, these definitions were 
based on flat plate geometry and this is quite reasonable for pipes 
except those having unusually thick walls. That is, in these defin- 
itions, pipe diameter was assumed to be quite large comparted to nine 
wall thickness. In subroutine AVERJ this assumption is not made. The 
appropriate definitions , based on truly cylindrical geometry , are given 
be lav, followed by demonstrations that they reduce to the code defin- 
itions as pipe radius goes to infinity while maintaining constant wall 
thickness. 

TAV, also designated as T, is the average metal temperature de- 



fined as 



as 


b 


r b 1 


T = /TdA//dA = 2 it 


Trdr/27T 


rdr = (1/hc) 


4 


a 


a 



b 



Trdr 



( 20 ) 



where a = inside radius , b = outside radius , h = b-a = wall thickness 
(do not confuse with earlier use of h as heat transfer coefficient), 



and c = (b+a)/2 = mean radius. If we define y = r-c, then 

f+h/2 r+h/2 

T = (1/h) (l+y/c)Tdy -* (1/h) 

J -h/2 •'-h/2 

and this is the definition given in the codes. 



rdy as c -*■ 00 



( 21 ) 
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We next define a quantity which can be termed the thermal 



+h/2 



moment , viz . : 

M = [r(T-T)dA = 2tt f r*(T-T)dr = 2ttc 2 [ ( 1+y/c) 2 (T-T)dy 

J J a ' -h/2 

Using the definition of T, this can be simplified to the form 



M = 2ttc 



f+h/2 

( 

-h/2 



( 1+y/c )iydy - 7rTh 3 /6 



( 22 ) 



(23) 



Tne next step in this development is to consider the average and 
the thermal moment of a linearly varying temperature distribution 
which has variation V from r = a to r = b. Such a temperature is 
given by T = Vr/h for which we find 



T = (V/h 2 c) 



rb 

r 2 dr = 



(V/3h 2 c)(b 3 -a 3 ) 



(PJD 



f h 

M = 2tt I r 2 ( Vr/n-T) dr = ( irVch 2 /6 ) ( 1-h 2 /12c 2 ) ( 2 r ; ) 

' a 

It requires considerable manipulation to obtain the last form. 

By equating these two values of M, one evaluates the (equiva- 
lent) linear variation 



V = [12 



C+h/2 

( 

-h/2 



( l+y/c)'.fydy - Th 3 /c]/h 2 (l-h 2 /12c 2 ) 



f+h/2 

(12Ai 2 ) Ttydy (26) 



J -h/2 

as c -*• The limiting values, as c -*• 00 , are, in each case, the 
definitions given by the various piping and vessel codes. The present 
subroutine, AVERJ, employs the more general definitions based on 
cylindrical geometry with finite radii. The variation V is also 
denotes as AT r in the codes and as CT1 in the output of the present 
program. 
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In the codes the quantity AT 2 is defined as the nonlinear 
residual temperature at the inner or the outer radius (whichever 
gives the greater result) after the constant part, T, and the linearly 
varying part, V, are subtracted. In the outout of the present program 
this quantity is denoted as DT2 and defined as 

DP2 = Mud | T^-T-DTl/2 1 , | T^T+ETl/2 1 } ( 27 ) 

The quantities AT t and AT 2 are intended, in the codes, to permit 
estimating significant thermal stresses in general pressure vessel 
components. However, in the particular case of uniform hollow cylin- 
ders , the stresses resulting from radial temperature gradient are 
definitely classified as "local thermal stresses" and there is no 
actual need to separate the actual temperature distribution as has 
been indicated above. Instead, the thermal stresses can be calculated 
directly, and this is done in subroutine AVERT. The quantities SL 
and SR are the values of circumferential (or axial — they are enual) 
stress due to the temperature distribution at inside and outside metal 
surfaces, respectively, under the assumption of zero net axial force 
and bending and twisting moment. The latter contributions may be de- 
termined separately by use of a conventional piping flexibility analysis. 

Because of the variability (with temperature) o^ the Young’s 
modulus and of the coefficient of thermal expansion, the calculations 
are not trivial. Discussion of this matter is given in appendix B 
hereto. 
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APPENDIX A Discussion of the subroutines 



Subroutines CSUBP, PRAND, and ADS MU receive as .Vtauts the tem- 
perature (°F) and the pressure (psia) of the steam flowing in the pine 
and produce as outputs the values of c p (specific heat at constant pres- 
sure), IJp (Prandtl's number), and p (absolute viscosity) o^ the steam. 

The calculations are by polynomials in T and P the coefficients of which 
were obtained by least squares analysis of data points read from the 
graphs on pages 293, 29^, and 297 of the ASHE Steam Tables ( ). It is 

implicitely assumed that the steam is saturated or superheated. If the 
steam is wet, the proper values are probably not obtained. 

Sixty number pairs were used in dtermining the coefficients in 
CSIBP, fifty-four pairs for PRAND, and twenty- four pairs for ABSMU. The 
only type of accuracy evaluation performed was to observe that each 
the 'subroutines was able to generate the innut number pairs within one 
or two percent (generally much closer than this). Such accuracy is 
greater than that implicit in convective heat tranfer theory and certainly 
at least as great as the accuracy with which the innut values were read 
from the graphs. 

Subroutine AITCH determines the film heat transfer coefficient h 
by use of the Colburn formula ( cf . Giedt ( 8 ) ) 

h = N ST Gc p (A-l) 

where G is mass per unit time per unit area, i.e.. 



G 



= n/A 



flow 



( A-2) 



and N om is the Stanton number, 

ui 
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(A-3) 



N = 0 023 M — 

ST PR fJ RE 

In this formula N PR is the Prandtl number and is the Reynolds number 

N FE = Gd/M (A->0 

where d is the pipe inside diameter. The Colburn fomula is one which 

is recommended for cases of well develooed turbulent flow at fairly 

hipti Reynolds numbers. It seems to be an appropriate choice except 

perhaps for cases in which the operating transient involves a very 

great decrease in mass flow rate. In such cases Q indeed approaches zero 

but not necessarily in the way represented by the Colburn formula. 

Incidentally, subroutine AITCH avoids a computational difficulty as 

G goes to zero by combining the equations in the form 

h = 0.023 C p G°* 8 (u/d)°* 2 / N 2 ^ 3 (A-5) 

Ihe specific heat c is evaluated at the fluid bulk temperature 

P 

but Npp and p are evaluated at the "film” temperature which is the 
arithmetic average of fluid bulk temperature and inside wall metal 
temperature. Since the latter is an unknown in the calculations, an 
iterative evaluation must be used. 

Subroutines METAL, COEXP, and EYUNG obtain pertinent properties 
of the metal pipe material corresponding to input values of the metal 
temperature. These subroutines presently provide data for only one 
particular material, the standard low chrome-mo ly high temperature 
alloy designated by the number 22 — pipe is designated as P22. Sub- 
routine METAL provides values of thermal conductivity (k), the rate 
of variation (k* = dk/dT), and the thermal diffusivitv using poly- 
nomial representations of material test data. 
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Although Table 1-4.0 of the ASME Nuclear Components Code (2) 
suggests using the same values for thermal conductivity and thermal 
diffusivity for P22 as for low carbon steel, data given in the TT’RC 
Data series, References (15) and (16) , indicate a sufficiently great 
dependence on alloying coastituents that we have preferred to use the 
latter. Thus we get substantially lower values of conductivity and 
diffusivity than the ASME tabulation gives. The actual data employed 
in constructing subroutine METAL are the data points from curve 19 
page 1158 of Ref. (19) and curve 2 page 343 of Ref. (16). It should 
be noted, however, that it would be a very simple matter to modify 
the subroutine METAL to use the ASME data rather than the TPRC data. 

Subroutines CQEXP and EYUNG are used (by subroutine AVERG) to 
provide the coefficient of thermal expansion and the Young's modulus 
corresponding to the given input temperature. The data is from 
Appendices D and C of the Power Piping Code, Ref. ( 3), and linear 
interpolation is used. Poisson's ratio is assumed to have the 
constant value 0.3. 

All subroutines as well as the main program use English custom- 
ary units, specifically inches (rather than feet) and seconds (rather 
than minutes or hours). 
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APPENDIX B Stress Calculations 



Assuming what is called a "quasi-static" situation, i.e., 
neglecting forces associated with accelerating the mass of the pipe 
as its thermal expansion causes displacement of its material particles, 
and assuming also that the difference between adiabatic and isothermal 
material constants is of neglegible significance, we nevertheless are 
faced with a nontrivial problem in determining the state of stress 
which results from the temperature changes discussed in the body of 
this report. 

Two material constants, E (Young’s modulus) and n(the coefficient 
of linear expansion — do not confuse with thermal diffusivity) are of 
significance in connection with stress evaluation. Both vary with 
temperature. Further, as will be discussed below, ot varies with stress 
state. A number of writers have dealt with stresses in tubes due to 
changes in temperature, among them Hilton ( 9) and Chang and Chu ( 7)* 



The analysis, in effect, reduces to the following. 

0 r = [E/r 2 ( 1-v ) ] [ ( r 2 -a 2 ) J(b)/(b 2 -a 2 ) - J(r)] (B-l) 

o - = [E/r 2 ( 1-v ) ] [ ( r 2 +a 2 ) J(b)/(b 2 -a 2 ) + J(r) - r 2 oT(r)] (B-2) 
0 

o = Ee + [E/(l-v)][2vJ(b)/(b 2 -a 2 ) - cYP(r)] (B-3) 

z z 

where r 

J(r) = f £«T(C)d£ (B-'l) 

‘ a 

Presuming that the axial force F is zero, we get 

g = 2 J(b)/(b 2 -a 2 ) (R-5) 

z 



and can thus evaluate a . For all temperature distributions which 

z 
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result from external heating (i.e., from heat applied to the inner or 
outer surfaces but not generated within the wall itself) the extremal 
stresses will be at r = a or at r = b. (With external insulation, it 
is plausible that the extreme must be at r = a.) We thus have = 0 
at r = a and at r = b, and we also have 

= [E/(l-v)][2J(b)/(b 2 -a 2 ) - ntT( a) ] at r = a (B-6) 

[E/(l-v)][2J(b)/(b 2 -a 2 ) - otT(b)] at r = b (B-7) 

In computing J(r) and J(b), the a under the sign of integration 
is the "instantaneous" value appropriate toT(r). It is (nresumably) 
the "zero-stress" value of a which is tabulated. Obviously, however, 
during the temperature changes the various elements of the nine wall 
are not stress-free and the question arises how it is that we 
should use the zero-stress value of a. The explanation is surely not 
well known (this question is not even considered in any of the standard 
references on thermal stress analysis); the explanation which follows 
was given by the writer in a report, dated November I960, entitled 
"Seflections on Piping Flexibility Analysis," which was intended to ex-^ 
plain the basis of procedures incorporated in the 'TOC Document" of 
the Pfechanical Design Committee of the ANSI (then ASA) B31 Piping Code 
Committee. 

This explanation considers a one dimensional constrained thermal 
expansion problem of a carbon steel bar of length L which at 70 17 
fits precisely between two rigid and immovable anchors; that is, at 70 F 
there is no stress and no gap. The temperature is increased to 300 F. 
One is to determine the resulting compression stress, oresuming that 
there is no tendency to buckle laterally. As fundamental we take the 
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relation 



de = do/E + otdT 



05-8) 



and, since, in this problem de = 0, we have 

da = -EadT (15-9) 



r300 



°F= - 



70 



E(T)a(T)dT 



where subscript („) indicates the "final" condition. 



( 55 - 10 ) 



Appendices Q and C of ANSI B31. 1-1973, Ref. ( ), give the 

following data for carbon steel pipe: 



T 


E*10 -6 


orlO 6 


Ea 


70 


27. 


6.01 


169.35 


200 


21.1 


6.38 


176.73 


300 


21. 1| 


6.60 


180. 84 



Using trapezoidal integration, we evaluate equation (B-10) to get 



° p = - ( 169 . 35+176 . 73 ) ( 130/2 )+ ( 176 . 73+180 . 84) ( 100/2 ) 



= -4037^ psi (B-ll) 

(This result is incorrect as we will show.) 

An alternate calculation obtains the answer by considering one anchor 
to be removed and letting the bar expand freely (asr*-s tress condition) 
from 70 F to 300 F and then supplying sufficient compressive load to 
restore the bar to its original length. This calculation gives 

„ f300 r 

°e = (T)dT = [(6. 07t6. 38)(65)+(6. 38+6. 60)(50)]-10 

> 70 

= 1.4583*10" 3 (B-12) 

a p = -(27. 4-10 6 )( 1.4583* 10~ 3 ) = -39956 psi (B-13) 

which is a slightly, but definitely, smaller value. 

Nov; it would indeed be disturbing if, without any inelastic 

behavior, the final stress state were to denend on the sequence of 
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load- temperature application. This would indeed be a new type of 
thermal ratcheting. However, we will new show that the second cal- 
culation gives the correct result for any load- temperature sequence, 
based upon the postulate that the final stress state is independednt 
of sequence and depends only on material nronerties, the final tem- 
perature state, and the final strain. In other words, we n resume that 
equation (B- B) reoresents an exact differential so that 

e p “ f(°p»Tp) 



If this is the case, we have 



9 / 3e _ 3 / 9e s 
lyracr “ 3cT3 r r 



and, since 



we have 



3e 

3T 



3e = 1 
3a E 



(B-liJ) 



(B-15a,b) 



3a _ d ( U 
3 a dTV 



(B-16) 

(v;e have written d rather than 3 on the right since, by its definition, 
E does not depend on a . ) 

Thus, if E varies with T, then a must vary with a; i.e., a = 
a(T.a). The tabulated values we use are values of a(T,0) and these 
were appropriate in the second numerical calculation but not in the 
first. Thus, for our original problem we should integrate 

da = -E(T)a(T,a)dT (B-17) 

taking into account the dependence of a upon a. We have 

ot(T,c) = a(T,0) + J^(|£)da = a(T,0) + a|*(|) (B-lB) 

so that 



da = -E(T)[a(T,0) + (a/E 2 ) (dE/dT) ]dT 
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= -E(T)ct(T,0)dT - adE/E 



(B-19) 



Thus , 



-a(T,0)dT = (Edo - odE)/E 2 = d(o/E) 



(3-20) 



Performing this integration, we get 




o p = (-lJ»583‘10" 3 )(27Jt-10 6 ) = -39956 psi 



as in the second calculation. 

In the more complicated, multi-axial stress situation with which 
this monograph is primarily concerned, the argument takes the following 
form. We imagine any element (dr,rd ,d%) of the pipe to be separated 
from the adjacent elements and to experience stress-free thermal ex- 
pansion resulting in extensional strains of the form 



Then surface forces are applied to the element to cause it to assume 
a form such that all elements can be fitted to gether in the final 
stressed state. Since equations (B-l) through (3-7) were derived by 
satisfying the equilibrium and constitutive laws, they indeed describe 
the final state. 

One final remark needs to be made concerning stress calculations 
for code purposes. Although the so called "simplified" methods of an- 
alysis provide for stress components corresponding to the temperature- 
like quantities AT, and AT 2 , it is very important to note that in the 
specific case of a uniform pipe or cylindrical vessel the stress state 
resulting from a radial temperature gradient does not have to be separ- 
ated into "secondary" and "local" parts but, instead, the actual stress 
may all be considered as "local." The word "local" is quite important 



e 




(B21) 
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in considering allowable states of stress. A definitive statement to 
this effect is to be found on the fifth page of the AS r E "Criteria" 
document. Ref. ( 4 ) , viz: 

"A special exception to these general rules is the case of the 
stress due to a radial temperature gradient in a cylindrical 
shell. It is specifically stated in 11-412 (m) (?) (6) of Section 
III, and in 4-112(1) (2) (6) of Appendix 4 of Division 2, that 
this stress may be considered a local thermal stress. In 
reality the linear portion of this gradient can cause defor- 
mation, but it was the opinion of the Special Committee that 
this exception could be safely made." 

Thus the stresses calculated according to equations (13-6) and (B-7) 
of this appendix themselves directly give the significant thermal 
stresses due to the radial gradient, and these are indeed local 
thermal stresses. 
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APPENDIX C Instructions , Data Preparation, Typical Problem 



Instructions for the use of the programs and preparation of in- 
put data are given, in the form of comments, in the initial portion 
of the program itself. This portion of the program reads as follows: 
Temperature distribution, etc., in uniform pipe with exterior 
insulation, convectively heated by fluid contents. The basic 
procedure is a finte difference treatment of the one (space) 
dimensional partial differential equation in cylindrical co- 
ordinates. Metal properties are assumed to vary with temper- 
ature. Surface heat transfer coefficient is determined by 
Colburn's relationship. The fluid content in this version is 
assumed to be superheated steam and the pipe material is P22. 

Input is by means of two kinds of cards. The first card reads 
OD, NT, TINIT, N, NIT, and NPROB according to FORMAT(3FiO.O ,315) . 
OD is the outside diameter in inches. NT is the wall thickness 
in inches. TINIT is the uniform initial metal temperature in 
degrees Fahrenheit. (For nonunifom initial temperature, see 
below.) N is the (even) nurrber of equal subregions into which 
the wall thickness is divided. (If input N is odd, it is round- 
ed up to the next even number.) NIT is the nurrber of iterations 
employed each time step in order to solve the nonlinear matrix 
equation. NPROB is an integer identification number. 

If TINIT < -500., the program looks for nonuniform initial metal 
temperature, reading enough cards according to FORMAT ( 7^10 . 0 ) 
to provide the (N+l) required initial temperatures. 
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Subsequent cards read I, NTEES, THU , TIM2,PHI1, PHI2, PRES1, 
PRESS, FL0W1, and FLDV/2 according to F0RMAT(2I5,8F8.0) . I is 
an identifying serial number which is useful for ordering the 
data but which is not used in the nrogram. NTEES is the nurrber 
of time subintervals into which the large time interval from 
time = TIKI until time = TIRE is to be divided. PHI1 and PHI2 
are values of the fluid temperature at times TIM1 and TIMS 
respectively, This tenperature is assumed to vary linearly 
with time. Similarly for PRES (fluid pressure in psig) and 
FIDW (fluid flow rate in pounds per second). 

Ihe data deck for a problem terminates with a card containing 
10000 in the first five positions. A nuntier of problems, each 
having a data deck as described above, may be stacked. The 
sequence of problems is normally terminated by a card contain- 
ing -1. in the first three positions. 

Some thought was given to providing an option which would permit 
several time intervals to be calculated per output line of print. 

It was decided that not only would this be a needless complication 
but also that frequent output is appropriate and useful during those 
periods when, for one reason or another, small time increments are 
employed. 

It is easily possible to make minor internal changes so that the 
nodal temperatures , rather than the results calculated in AVERJ, 
are printed. This is useful when comparing results with those ob- 
tained by other methods of calculation. 
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The program listing is not given here for reasons of space economy. 
Those who are seriously interested in installing the nrogram should 
communicate directly with the writer to make appropriate arrangements . 
The program is written in FORTRAN using none other than quite primi- 
tive capabilities. It has been run on an IBM 360/67 requiring only 
56 K core for execution (both G and 11 compilations). For a stack o r 
thirty two problems, involving a total of 9240 time steps, with 1 = 

10 (10 layers) and NIT = 5 (five iterations), time requirements were 



as follows: 




G Compiler 


H Compiler 




Compile 


14. 38s 


23.20s 




Link-edit 


1.87s 


1.74s 




Execute 


I4ml9.48s 


llml7.03s 




Total 


1%i35.73s 


llm4 1.97s 
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Figure C-l Instructions for the use of PIPETEM 
This material is a oart of the PIPETEf'’ nrorrram itself 
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APPENDIX D Program verification 



By artificially causing the program PIPETEM to operate with constant 
values of diffusivity, conductivity, and surface heat transfer coefficient, 
it may be used to solve problems for which "classical" solutions are known. 
The program was initially verified by taking O.D. = 10000., W.T. = 1., and 
comparing results with data points taken from graphs given by McNeill and 
Brock (12). When a number of such comparisons checked within the accuracy 
with which the graphs could be read, a more exacting test was made. This 
is described in what follows. 

Ozisik (13) gives solutions involving true cylindrical geometry. v or a 
uniform cylinder with inside radius a and outside radius b, insulated on 
its outer surface, and convectively exposed, with constant surface heat 
transfer coefficient h, to an internal fluid having "ramp" temperature 
= At (where A is a constant and t denotes time) , and having initial 
temperature distribution T(r,0) = 0 at time t = 0, the solution may be 
written as 

TCr.to - l x m < t)* (p)/y (l) 

m=l 

where 

x (t) = 2ABot[exp(-y t) - 1 + y t]/a 2 
m m m 

Y m = “ ( V ' a)2 

yo> - voy/ysy - j .(oy/ j ,<oy 

= ytc^ysj/yi)] 1 - i -<B/yn 

B = ahA (Biot's number), p = r/a, o = b/a 

and the numbers y are the zeros of the function 
m 

[BJ 0 (y) + yJ 1 (u)]Y 1 (py) - [BY 0 (y) + yYj (y) ]Jj ( oy) 

The symbols J and denote Bessel functions. 
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It proved to be quite difficult to devise a program to nerform 
these calculations. An apparently different solution, given by Luikov 
(10), was dealt with first, but debugging was not completed successful!/. 
A later work by Luikov (11), in Russian, nave the form we actually used, 
citing ftzipik (13) as source. Both single and double precision versions 
of an implementing program were constructed and debugged. The double 
precision version, much slower than its counterpart, partly because of 
the use of a "home-made" DP Bessel function subroutine, was used to 
verify that the single precision version provided adequate accuracy. 

In all such checks, results from the two programs agreed to five sig- 
nificant digits. 

The test problem employed to check the integrity and accuracy 
of PIPETEM was as follows : a = inches, b = 5 inches, A = 1. , B = 1. , 
a = 1., k = 1. , k* = 0., h = 0.25. Analytical calculations were made 
for t = 0, 1., 2., 3.> and 5. PIPETEM calculations were made for 
t = o. to t = 5. at intervals of 0.1. In PIPETEM we took N = 10 and 
HIT = 5. (V/e wasted a few seconds by not usinf HIT = 1. This problem 
is completely linear and no iteration is necessary.) In the analytic 
solution initially we took account of ten terms in the surrmation. 

The results were so close that differences can not be displayed 
graphically. However, the differences were sufficiently great as to 
demand an explanation. Truncation error in the analytic solution was 
suspected. Accordingly, the analtic solution was repeated three more 
times using five, twenty, and forty terms. (The forty term solution 
was very extravagant of computer time.) Plotting typical analytic 



results versus reciprocal number of terms clearly indicated that Indeed 
the analytical solution did suffer from truncation error and that extra- 
polation did indeed lead to the result given by PIPETEM. Bince the theo- 
retically "exact" solution has by far the greater error, it is not easy 
to estimate the accuracy of the PIPETEM results, but in this problem the 
error appears to be only littl« n«reth?n one part in one thousand. The 
following tabulation gives calculated results for T(0,5), for example. 

N ■ 5 N = 10 N = 20 N = 40 N = 00 PIPETEM 
2.05694 2.08627 2.09991 2.10647 2.11496 2.11220 

Table D -1 Calculated results for T(0,5) 

The first four tabulated results were obtained by the single precision 
analytic solution. (The double precision version also gave 2.0R627 for 
N =10.) The extrapolated value, for N = 00 , was obtained by assuming 
that T = A + Bx + Cx 2 + Dx 3 , where x = 5/M, leading to the evaluations: 

A = 2.114958096, B = -.073526667, C = .049093333, D = -.033584762. The 
difference between the extrapolated analytic value and the value provided 
by PIPETEM is 0.00275 which is 0.13% of the average of these values. 

As an ultimate check on the validity of the program, an artificial 
nonlinear problem was made up and solved numerically. The conditions 
are: inside radius a = 1, outside radius b = 2, conductivity k = T (so 
that k' =1), heat transfer coefficient h = 15e^, fluid temperature <t> 

= 6e^, initial temperature distribution T(r,0) = r + 4/r. The diffus- 
ivity a is a function of both temperature and time, viz: 
a = Y/[l+( 4-Y) 2 /(4+Y) 2 ] , Y = (Z-/Z 2 -l6) 2 /4, Z = Te -t 
It is easy to verify that T(r,t) = e fc (r + 4/r) is a solution of this 
problem. Althoupfc matters of uniqueness of solution to nonlinear sys- 
tems are difficult to deal with, it happens that PIPETEM indeed yields 



the solution cited. 
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Attention was focussed on the values = T(l,t*) = 100.000 
Tg = T(2,t*) — 80.000, where t* = log^20. Percent errors are 
shown in Table D-2 corresponding to the 128 evaluations indicated, 
values of r being 1, 2; values of II being 2, 4, 6, 10; values of 
NIT being 1, 2, 5, 10; and values of DTFE being t*/5, t*/20, t*/50, 
t*/200. 



DIEE NIT N = 2 



!I = 4 



II = 6 



N = 10 



t*/5 



t*/20 



t*/50 



t*/200 



1 

2 

5 

10 

1 

2 

5 

10 

1 

2 

5 

10 

1 

2 

5 

10 



-21.5V-l4.04 
-2.726/+1.479 
-2.592/-6.721 
-.8019/+4.055 
-2.042/+2.726 
-1.962/-. 9932 
-1. 648/+. 8822 
-1.735/-. 1605 
-1.379/+2.520 
-1.840/-. 1536 
-1.714/+.4576 
-1. 723/+. 3710 
-1.631/+1.020 
-1.73V+.3131 
-1. 725/+. 3715 
-1. 725/+. 3713 



-14.70/-7.173 
+.3159/+3.258 
-l.78V-6.35l 
+.5816/+3.232 
+1.085/+4.366 
-1.070/-1.907 
-.3848/+. 7397 
-.5519/-. 0690 
+.2103/+2.209 
-.7005/-. 556 8 
-.5162/+. 1779 
-.5285/+. 0778 
-.3372/+. 7247 
-.5 426/+. 0234 
-.5287/+. 0857 
-.5288/+. 0489 



-9.532/-2.487 
+.0102/-. 0009 
-.9409/-3.697 
-.0216/+. 3814 
+1.840/+4.425 
-1.010/-2.239 
-.1063/+. 5389 
-.2782/-. 0115 
+.6620/+2.185 
-.4543/-. 5960 
-.2334/+. 1084 
-.2466/-. 0069 
-.0229/+. 6389 
-. 2618/-. 0232 
-.2464/+. 0365 
-.2465/+. 0351 



-3.372/+3.670 
-.7222/-2.646 
+1.289/+3.447 
- 2 . 023 /- 6.768 
+2.569/+4.810 
-l.10V-2.594 
+. 067 1/+. 5541 
-.1289/-. 1902 
+.9698/+2.209 
-.3513/-. 6870 
-.0781/+. 0837 
-.0916/-. 0043 
+.I658/+. 6 182 
-.1089/-. 0567 
-.0910/+. 0131 
-.0911/+. 0114 



Table D-2 Percent errors in evaluation of T^ and Tg in test problem. 
Error in T A (Tg) appears to left (right) of slash (/). 



It may be noted that the accuracy is quite sufficient for engineering 
purposes if the number of layers N = 4, the number of iterations NIT 
= 5, and the total time interval is divided into 20 parts. It may 
also be remarked that the total time to evaluate all 128 results was 
only 2 minutes 20 seconds. 



- 27 - 



APPENDIX E Heat Transfer at Outer Radius 



Following the completion of the preceding portions of this renort, 
an application came to attention which called for the aumentation of 
PIPETEM capabilities. Ordinarily high temperature pining and fittings 
are lagged (insulated). However, in some cases it might be of inte v ’- 
est to follow the temperature history of unlamged pipes when subjected 
to internal fluid transient temperatures. Also, PIPETEM is obviously 
applicable to insulated flanges if one discounts axial conduction and 
the additional thermal resistance of the flange-bolt interfaces. 
Down-transients (i.e., decreases in the temperature of fluid contents) 
lead not only to local thermal stresses but also to discontinuity 
thermal stresses if there are significant changes in wall thickness, 
since the average temperature of the thicker part lags (in tine) be- 
hind that of the thinner part. This is particularly tree of flange- 
pipe assemblies. Leaving the flange unlagged (but with lagging on 
the pipe itself) causes the average temperature of the flange to be 
lower than that of the pipe since the former can dissipate heat by 
convection and radiation. Then, when a down- transient occurs, the 
resulting thermal dislocation stress will be less than what it would 
have been if the flange had been lagged. The calculation of the 
temperatures in an unlagged flange really calls for finite element 
treatment since axial conduction and heat loss r rom bolts, nuts, and 
flange face is obviously of importance. However, a useful estimate 
of the effect can be gained by considering only radial flow, with, 
perhaps, an artificially enhanced value of the effective convective 
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external heat transfer coefficient. 



Accordingly, PIPETEM has been modified to account for a convective 
condition at the outer boundary, viz. 

k^3 = hOHJ? M ) E-l 

rather than B = 0. The ilth finite difference equation becomes 

t6 ri + 2 + a+ C2 + i - V + (< r » + ) )><, - rf, 

= <S m -2)T h + 2T n - o(T k -*)[2 + i - nn(- M -',) J 

+ cV(2 + i + pJpV) E-2 

This equation, which takes the place of equation 16 and which reduces 
to equation 16 when h (external) = 0, my be compared term by term with 
equation 14. In equation E-2, quantities have the meanings ascribed 
to them on page 5, but the evaluations are at the outer radius. Note 
that all signs are the same (as in equation l l \) excent for the i terms. 
The quantity is the external ambient temperature (for convection). 

External heat transfer can take place both by convection and by 
radiation, so that the coeffcient h in equation E-l is the sum of 
two coefficients 

h = h c + h r E-3 

The first, h , is the actual convective coefficient. We have presumed 
c 

natural convection from a cylinder into still air, using the relation 

h = 0.2>l(e/d)°* 2Ij E-'< 

on page 218 of reference 8, Here 

0 = |T m - E-5 

and d is the diameter of the cylinder, i.e., the outer diameter of the 
pipe. The "radiative-convective" coefficient h r is an equivalent con- 
vective coefficient to account for radiative heat transfer and is de- 
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fined by the equation 

Q = aeA [ ( Tpij+C ) (x^J+C) ** ] = h^Tg-iJO E-6 

in which a is the Stefan-Boltzmann constant (0.171A xlO -8 Btu/hr ft 2 P 1 *), 

€ is the emissivity (we take 0.8 for bare steel), A is area (which cancels), 
iJj is the temperature of the surrounding with which radiative heat trans- 
fer takes place (this may be and usually is different from the convective 
ambient temperature 40, and C = d6o is the constant required to convert 
Fahrenheit temperatures to Rankine (absolute) temperatures. 

The value of h determined from equation E-3 is in units of feet and 
hours rather than inches and seconds and it must be multiplied by 1/(720) 2 
before being used in PIPETEM. There is provision in PIPETEM for using 
the relation 



h = f h + f h 
c c r r 



E-7 



where f and f are "augmenting factors" which may be used to account 
for the fact that the geometry is not exactly that of an infinite cylinder 
so that a "face" or protuberances such as nuts and bolts rnav convect 
and/or radiate more heat than would otherwise be accounted for. 

The convective ambient temperature ip and the radiative ambient 
temperature i p could be made to be functions of time, ,1ust as with the 
internal fluid temperature <t>, but this has not been done in the present 
inplementation in which and ip are taken to be constant. 

The input (and the input instructions) have been modified slightly 
to accorrmodate this added facility. For each problem an additional card 
is inserted after the first card described previously. If this card is 
blank, perfect insulation is assumed at the outer surface. Otherwise , 
the card inputs \p t ip, and c. 
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Another slight modification has also been made in PIPhTEH. In 
the first card of a problem deck, settinm TINIT = - r 500. causes the 
program to search for and read initial temperature data as described 
earlier, while setting TINIT to any value greater than -485. causes 
the program to assume a constant initial tennerature enual to the 
value of TINIT. The new capability is that setting TINIT = -400. 
causes the program to start with initial temneratures enual to the 
final temperatures for the preceding problem. Unless the values of 
N are equal for the two problems, some nonsense is generated. In 
use, of course, this capability is intended for a succession of two 
or more problems involving the same pipe, divided into the same num- 
ber of layers. A modified set of COMMENT cards at the beginning of 
the program describes the modified capabilities and the details of 
the problem deck. 
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